Computer Simulation Study of the Phase Behavior and Structural 
Relaxation in a Gel-Former Modeled by Three Body Interactions 

Shibu Saw 1 , Niels L. Ellegaard 1 , Walter Kob 2 , Srikanth Sastry 1 
1 Theoretical Sciences Unit, Jawaharlal Nehru Centre for Advanced Scientific Research, 

Jakkur Campus, Bangalore 560 064, India. 
2 Laboratoire des Colloides, Verres et Nanomateriaux, 
UMR5587 CNRS, Universite Montpellier 2, 34095 Montpellier, France 

(Dated: January 20, 2011) 

Abstract 

We report a computer simulation study of a model gel-former obtained by modifying the 
three-body interactions of the Stillinger- Weber potential for silicon. This modification re- 
duces the average coordination number and consequently shifts the liquid-gas phase coex- 
istence curve to low densities, thus facilitating the formation of gels without phase sepa- 
ration. At low temperatures and densities, the structure of the system is characterized by 
the presence of long linear chains interconnected by a small number of three coordinated 
junctions at random locations. At small wave- vectors the static structure factor shows a 
non-monotonic dependence on temperature, a behavior which is due to the competition 
between the percolation transition of the particles and the stiffening of the formed chains. 
We compare in detail the relaxation dynamics of the system as obtained from molecular 
dynamics with the one obtained from Monte Carlo dynamics. We find that the bond cor- 
relation function displays stretched exponential behavior at moderately low temperatures 
and densities, but exponential relaxation at low temperatures. The bond lifetime shows an 
Arrhenius behavior, independent of the microscopic dynamics. For the molecular dynamics 
at low temperatures, the mean squared displacement and the (coherent and incoherent) 
intermediate scattering function display at intermediate times a dynamics with ballistic 
character and we show that this leads to compressed exponential relaxation. For the Monte 
Carlo dynamics we find always an exponential or stretched exponential relaxation. Thus 
we conclude that the compressed exponential relaxation observed in experiments is due to 
the out-of-equilibrium dynamics. 
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I. INTRODUCTION 



Gels are ubiquitous in daily life including food, cosmetics and medicines. They are low 
density disordered networks of interacting molecules which can sustain weak stress. In this 
sense they behave like solids^. In spite of their low density, fluids that form gels exhibit 
slow relaxation dynamics at low temperatures, similar to glasses. Based on the life time of 
bonds between constituent units, gels can be classified as chemical or physical gels^-. While 
chemical gels are formed by the formation of strong covalent bonds in a disordered structure, 
physical gels (such as gelatin, colloidal gels, etc.) arise due to relatively weak interactions. 
The mechanisms for the formation of such physical gels, i.e. gelation, are still not fully 
understood. Gelation can arise in phase separated solutions^ 1 ^ due to the intersection of 
the glass transition and spinodal lines^, by the formation of a system spanning network 
in homogeneous suspensions^ 1 ^, or due to aggregation of small clusters by the process of 
diffusion limited cluster aggregation^^. 

Gels exhibit unusual dynamical properties. Their relaxation process is often stretched 
exponential^, a compressed exponential^ 1 ^-, or logarithmic^. Experiments on non- 
equilibrium colloidal gels show compressed exponential relaxation behavior of the dynamic 
structure factor with an exponent of ~ l.f>= How these features are related to the equilib- 
rium properties of the gel-former and to the details of the interactions between the particles 
is presently not known. Unfortunately there exist so far only relatively few studies of sus- 
pensions in equilibrium which approach gelation while remaining homogeneous^ 1 ^ 1 ^^— . 
In order to obtain such gels it is necessary for the system not to enter the liquid-gas (LG) 
coexistence region on its approach to the gelation line. This can be accomplished by shifting 
the liquid-gas coexistence region to lower densities and temperatures, thereby opening up a 
low temperature, low density region of the phase diagram where the solution remains homo- 
geneous but is characterized by the presence of bonds which are strong relative to thermal 
fluctuations. One way to effect such a shift is to reduce the maximum coordination number 
which a particle can have^. A shifting of the LG coexistence region in this manner has 
been realized for patchy colloidal particle suspensions (where colloidal particles have a small 
number of attractive patches on its surface)^2i, in maximum valency lattice gas models^, 
hard sphere models with short range square well attractive interactions^, and models with 
dipolar interactions^. Although this scenario has not been routinely realized experimentally, 
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a recent study of colloidal clay (laponite) presents an experimental realization of the above 
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scenario^. 

The reduction in maximum coordination number can also be obtained using three body 
interaction potentials^. Three body interactions are expressible in terms of angles formed 
by triplets of particles, and can be chosen to energetically stabilize open network structures 
in which the particles have very small (tunable) connectivity. One recent example of this 
approach has been presented ivM-, wherein the Stillinger- Weber (SW)^ model potential for 
silicon has been modified to generate structures with low coordination. 

In the present paper, we study the static and dynamic properties of gel-forming fluids 
obtained by the above-mentioned modification of the SW interaction potential. Note that 
all the presented results are for the system in equilibrium, i.e. they are not affected by aging 
phenomena. The remaining of the paper is organized as follows: In Sec. II, we describe the 
SW model potential of silicon and its modification for homogeneous gelation, and provide 
the computational details relevant to this work. In Sec. Ill, we discuss the static properties 
of the system and the dynamic properties in Sec. IV. Finally, Sec. V contains a summary 
and conclusions of our work. 

II. MODIFICATION OF STILLINGER- WEBER POTENTIAL AND COMPUTA- 
TIONAL DETAILS 

Here we describe the SW— potential which we modify in this work. The SW potential 
was originally proposed with parameters chosen to provide a reasonable description of the 
experimentally observed thermodynamic and structural properties of crystalline and liquid 
silicon. Its functional form is given by two-body and three-body interaction terms and is 
written as 

usw = ^2v 2 (r ij /a) + ^ ^(r^/V, r^/a, r fe /cr), (1) 

i<j i<j<k 

where a is the diameter of the particles, rj is the position of particle i, and is the distance 
between particles i and j. The two-body potential is short-ranged and has the form 
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, Ae(Br~ 4 - 1) exp (^-) r <a 

war = { > (2) 

> a 

where A = 7.049556277, B = 0.6022245584, and a = 1.8. The repulsive three-body potential 
is also short-ranged, and is given by 



v 3 (r u Tj, r k ) = h(rij, r ik , 9 jik ) + h{r ih r jk , 9 ijk ) 

+h(r ik ,r jk ,9 ikj ), (3) 

where 8ji k is the angle formed by the vectors and r ik and 

h(rij,r ik , Ojik) = e\ exp[ — 1 — — ] (cos 9 jik + a) 

Tij d T{ k (1 

xH(a-r ij )H(a-r ik ), (4) 

where A = 21.0,7 = 1-20, and H(x) is the Heaviside step function. The choice a = 1/3 in 
(cos 9 ji k + a) 2 favors a tetrahedral arrangement of atoms as found in silicon. 

While the two-body interaction favors a close packed arrangement of particles, the three- 
body interaction term favors an open structure whose geometry depends on the parameter a. 
The balance between the two-body and three-body interactions, controlled by the parameter 
A dictates the final preferred geometry of particle arrangements^. The idea behind the 
modification of the SW potential is to adjust the three-body interactions so that the average 
coordination number is reduced, which in turn will shrink the phase-coexistence curves 
and thus increase the region in the T — p diagram in which gels can form, without the 
involvement of phase-separation. Although similar in spirit to models, e.g. in&2i, in our 
model the number of bonds depends on density and temperature, which is more similar in 
this respect to the model studied irA 

By appropriately choosing the values of A and a, it is possible to obtain a system that 
has only a small coexistence region, with a network morphology whose average connectivity 
can be tuned. We choose (among other possible choices) values A = 10 and a = 1.49 since 
these values avoid also the formation of ordered structures (further details may be found 
in Ref2&22). In the following, all quantities are reported in reduced units for the modified 
Stillinger- Weber potential, i.e. e, a and m are the units of energy, length, and mass. 
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We have performed molecular dynamics (MD) simulations in the NVT ensemble using 
4000 particles in a cubic box with periodic boundary conditions. We have used the con- 
straint method^ for constant temperature simulation, where the equations of motion are 
modified to maintain a constant kinetic energy. The MD simulations for the SW poten- 
tial are computationally demanding due to presence of the three-body interaction energies, 
where angles of all the triplets need to be determined. Therefore, we have used an efficient 
method which allows the calculation of the three-body interaction term by means of a two 
loop summation^— . We have extended this approach to calculate the force due to three- 
body interactions, as we describe in the Appendix. We have used a time step of At = 0.005 
in the MD simulations. The temperatures investigated are T = 5.0, 3.0, 1.0, 0.30, 0.20, 
0.10, 0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.03, and 0.028 at density 0.06. Most of the data are 
averaged over 5 different samples. We have used 20 to 30 million MD steps for equilibration 
and 30 to 100 million MD steps for data generation. 

The Monte Carlo (MC) simulations have been performed with 512 particles, typically 
using a maximum step size of 0.085. In the MC simulations, 2 million MC steps were used 
for equilibration and 10 to 20 million MC steps for data generation. 

III. STATIC PROPERTIES 

Figure [T] shows the radial distribution function, g(r), for different temperatures at density 
0.06. At higher temperatures g(r) shows a single peak corresponding to the interparticle 
distance. Additional small peaks appear at lower temperatures corresponding to second and 
third nearest neighbor distances, similar to dense liquids. However, here the g(r) curves 
are characterized by relatively sharp peaks, with the gaps in between having values of g{r) 
close to the ideal gas value of 1 (except for low temperatures, between the first and second 
peaks of the g{r)). These features arise from the fact that the system is composed, at low 
temperatures, of long interconnected chains. 

In order to study the structure on large length scales, it is useful to consider the structure 
factor, S(k), defined as 

1 N N 

^ k ) = A7EE( ex p(- ik -^- r '))) ' ( 5 ) 

j=l 1=1 
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FIG. 1: The radial distribution function, g(r), at density p = 0.06 for different temperatures. At 
high temperatures, g(r) does not show much structure beyond the first peak. At low temperatures 
g(r) exhibits multiple peaks corresponding to well defined separation distances for second, third, 
... nearest neighbors. 

where N is the total number of particles and k is the wave- vector. Figure [2] shows the S(k) 
for a range of temperatures at density 0.06 and 0.10. 

At high temperatures, S(k — > 0) has a value close to the ideal gas value of 1, but as 
the temperature is lowered, the S(k) at small k values gets substantially larger—. Note 
that at intermediate and small wave vectors, S(k) shows a non-monotonic dependence on 
temperature, displaying a maximum value at intermediate temperatures (see inset in Fig. |2K 
where we show the T dependence of S(k m i n ), i.e. the value of S(k) at the smallest wave- 
vector k m i n that is compatible with the size of the simulation box). For p = 0.06, the 
peak occurs at T = 0.08, i.e. well below the percolation transition at T = 0.115^, where 
one may naively have expected the maximum to lie. The reason for the shift away from the 
percolation point may be rationalized by noting that in addition to percolation, which would 
give rise to a power-law in S(k), the structure of the system changes significantly because 
of the formation of increasingly longer and stiffer linear chains. The latter effect will, via 
the form factor of the chains, also give rise to an increase in S(k) at low k. The net effect is 
a shift of the temperature at which S(k) has the maximum value at low wave vectors (see 
phase-diagram in Ref.— ). We have calculated S(k) for system sizes N = 512, 4000, and 8000 
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FIG. 2: The structure factor for different temperatures at (a) p = 0.06 and (b) p = 0.10. The S(k) 
at the lowest wave-vector k shows non-monotonic behavior in T for both densities, as shown in the 
inset of panel (a). Temperatures below the maximum of S(k m i n ) are shown in solid symbols while 
those above are shown in open symbols to clearly indicate the non-monotonic behavior of S(k). 

for some state points in order to investigate finite size effects and do not find any significant 
size effects in the form of S(k). 

The density dependence of S(k) at T = 0.04 is shown in Fig. El At low densities, 
p < 0.12, the first peak is at k — 6.15 corresponding to the nearest neighbor distance. 
With an increase of the density, a peak around k = 3.0 develops, becoming pronounced at 
the highest density shown; the value of S(k) at the smallest wave-vectors gets suppressed, 
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FIG. 3: The density dependence of S(k) for T = 0.04. For densities above 0.1, the large values of 
S(k) at low k get suppressed, and a finite k peak develops around k = 3.0. 

indicating a reduction in the compressibility and density fluctuations on long wavelengths. 

The local geometry of the system can be characterized by the distribution of the coor- 
dination number n. Two particles are considered to be neighbors if their distance is less 
than the (T and p-dependent) location of the first minimum of g(r). Figure H] shows the 
temperature dependence of P(n), the fraction of particles having coordination number n, 
vs. 1/T for density p = 0.06 (filled symbols) and p = 0.10 (open symbols). For n = 
and 1, the P(n) decreases rapidly and monotonically as T decreases, whereas for n = 2 the 
P(n) increases monotonically with decreasing T. These changes correspond to the forma- 
tion of a network of particles at the expense of isolated particles or dimers. Interestingly 
for n — 3 and 4 the T— dependence of P(n) is non-monotonic, displaying a maximum at the 
temperature at which S(k m i n ) displays a maximum (see Fig. |2]). Above this temperature 
of maximum P(3) and P(4), the structural change upon lowering T arises from a growth 
in the number of bonds leading to percolation. In contrast, at temperatures below the lo- 
cation of the maximum, the system forms increasingly longer linear chains at the expense 
of cross-links with larger coordination. In contrast, in the case of the gel-forming model 
irA P(3) increases with lowering of T and saturates at very low temperatures. The fraction 
of three-fold coordinated particles, P{3), increases strongly with density (data not shown). 
Thus at higher density there are more anchor particles and consequently the lengths of the 
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linear chain segments are smaller compared to that at low density. 

We have also calculated the cluster size distribution for density 0.06^. Even for tempera- 
tures near the percolation transition, we do not find a distinguishable power-law regime in the 
cluster size distribution, in contrast to the findings for other gel forming systems for which 
a power-law dependence corresponding to random percolation has been observecr^i'^. We 
believe this to arise from the fact that while the (exponential) cluster size distribution is ob- 
tained from a random aggregation process at high temperatures, the percolation transition 
is influenced by the emergence of increasingly long linear chains. 

The distribution of segment lengths shows an exponential behavior (see Fig. E]), with an 
average segment length that increases with decreasing temperature, as can also be inferred 
from the T-dependence of P(2) and P(3) in Fig. HJ In the inset of Fig. Owe show the average 
segment length obtained from an exponential fit to the distribution of segment lengths, as 
well as the ratio P(2)/P(3), which provides an estimate of the average segment length. 
These two are in good agreement. For other densities also, the segment length distribution 
exhibits exponential behavior and the average segment length decreases with an increase of 
density. 




-| Q I i I i I i I i I i I i I i I i I i I i 

4 8 12 16 20 24 28 32 36 40 
1/T 

FIG. 4: The T-dependence of P(n), the fraction of particles having coordination number n, at 
density 0.06 (solid symbols) and 0.10 (open symbols). At low temperatures most of the particles 
have coordination n = 2. With the decrease of temperature all coordinations decrease except the 
one for n = 2. 
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FIG. 5: The segment length distribution at density 0.06 for different temperatures, which exhibit 
exponential behavior. The lines are exponential fits to the data. Inset: The average segment length 
vs. 1/T from exponential fits in the main panel, and the estimate using the ratio of number of two 
and three coordinated particles, P(2)/P(3). 

The squared end-to-end displacement (S 2 (N C )) of chain segments shows a quadratic be- 
havior for small chain lengths N c and is linear for large N^. The persistence length can 
be denned as the maximum length of a chain for which (5 2 (N C )) ~ N% holds. At T = 0.028 
and density 0.06, the persistence length is 15 and it decreases with increasing temperature 
or density. 

IV. DYNAMICS 

In order to study the dynamics of the system we first consider the mean squared dis- 
placement (MSD) of particles, defined as 



which is shown in Fig. [(2 for different T at density 0.06. 

At high temperatures, a smooth crossover is seen between the short time ballistic and 
long time diffusive regimes. If the temperature is lowered, one sees the emergence of incipient 
plateaus, as seen in dense supercooled liquids. However, unlike dense liquids, one observes 
the presence of two "shoulders", rather than one^. The first appears at a length of ~ 0.4<r, 
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FIG. 6: Mean squared displacements for different temperatures at p = 0.06 in a log-log plot. The 
location of the I s * and 2 nd shoulders, seen for low T and discussed in the text, are marked. The bold 
dashed line shows the MSD from constrained simulations at T = 0.04 (see text) which saturates 
at the second shoulder. The inset shows the MSD at the lowest T for different densities. 

analogous to what is seen in dense liquids. The second shoulder appears at ~ 7a, indicating 
that transient localization occurs here at a much larger length scale compared to dense 
liquids. Although large, this length is significantly smaller than the average segment length 
of the chains (around 100 at the lowest T; see inset of Fig. [5]). Thus, we conclude that 
the amplitude of floppy motions contributing to the MSD at the second shoulder is much 
smaller than the length of the chains. 

The inset in Fig. [6] shows the evolution of the MSD at T = 0.04 with density. As 
density increases, the displacement corresponding to the second shoulder decreases and also 
it becomes less pronounced since the segment length between anchor particles decreases. 

Also included in the main panel of Fig. [6] is the MSD as obtained from constrained MD 
simulations at T = 0.04, i.e. a dynamics in which bonds are prevented from breaking or 
forming^. As we see in the figure, the MSD from constrained MD simulations saturates at 
the second shoulder, indicating that the caging associated with the second shoulder is that 
of chain segments executing floppy motion between two anchor particles, whose magnitude 
is constrained by the network topology. 

The non-Gaussian parameter has been studied extensively^ in the context of supercooled 
liquids as a way of characterizing dynamical heterogeneities. The non-Gaussian parameter 
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at2(t) is defined as^^- 



a 2 {t) 



3(r 4 (t)) 
5(r 2 (t)y 




(7) 



where 



( r4 w> = af£ k*)- r <(°) ) 



(8) 



t=l 



Figure [7Ji shows a^t) for different T at density 0.06. At high temperatures, a single peak 
in a^{t) is observed at around t — 10, a time that corresponds to the crossover in the MSD 
from ballistic to diffusive behavior. At intermediate temperatures, an additional peak at 
t « 1 emerges, which is the time at which the first shoulder in the MSD is observed. Its 
origin is likely the heterogeneous dynamics related to the fact that different particles can have 
quite different type of cages since they can be two or three-fold coordinated. If T is decreased 
even further, the location of the peak at longer times rapidly shifts to larger times, and its 
height increases with decreasing T. This peak corresponds to the second shoulder seen in 
the MSD. We note that the height of this peak is, even at the lowest temperatures, less than 
0.2, a value that is significantly smaller than the ones that are found in dense glass-forming 
systems, such as, e.g. Lennard- Jones mixtures^ for which a2{t ma x) ~ 1-5. This difference 
is probably due to the fact that in the present system there is, at low temperatures, not 
that much variance in the relaxation dynamics of the individual particles since the main 
relaxation process is a cutting of a chain and a subsequent reconnection of the loose ends to 
the rest of the network. This process depends only very weakly on the particle considered 
and hence the ^(i) is small. 

In Fig. [7b we show the density dependence of the non-Gaussian parameter for T = 0.04. 
The position of the first peak, around t — 1, is independent of density whereas the second 
peak position moves to larger times with decreasing density, in agreement with the density 
dependence of the MSD on the time scale of the second shoulder. 

To study the life time of the bonds, we calculate the bond correlation function (psif) 



defined as^ — 




(9) 



where 
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FIG. 7: The non-Gaussian parameter a% (t) for (a) different T at density 0.06 and (b) different 
densities at T = 0.04. 

Sriijit) = nij(t) - (n) 

and riij = 1 - H(r i:j - r cut ) . (10) 

H(x) is the Heaviside function and r cut is the bond length defined as the distance at which the 
first minimum of g(r) occurs, riy is 1 when a bond is present and otherwise. Hence, </>s(t) 
counts the fraction of bonds found at t = that survive after a time t, independent of any 
breaking and reforming at intervening times. Figure |8] shows the bond correlation function at 
density 0.06 for different temperatures (upper curves) and at T = 0.04 for different densities 
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(lower curves). (For the sake of clarity, the lower curves have been shifted down by a factor 
of 100.) In order to recognize better the T and p— dependence of this correlation function 
we plot it as a function of t/rs, where t# is the average bond life time as determined from 
the area under (f>B (t). We see that the bond correlation functions decay exponentially with 
time at low T, but that the decay is slower than exponential at intermediate temperatures. 
At high T the decay is again exponential (data not shown). Similarly, at T = 0.04, the 
decay behavior is exponential at low densities and stretched exponential at higher densities. 
In both cases, the stretched exponential behavior appears to be associated with an increase 
of disorder in the local environment of the bonds. 




FIG. 8: The bond correlation function (j>B(t) as a function of t/rs at different temperatures and 
density 0.06 (upper curves ) and at different densities for T = 0.04 (lower curves). The latter set 
of curves has been shifted down by a factor of 100 for the sake of clarity, tb is the average bond 
life time. 

Figure [9] shows the T-dependence of the bond life time, tb, for density 0.06 and 0.10 from 
MD simulations as well as the tb from MC simulations at density 0.06. We recognize from 
the figure that the bond life times from MD and MC simulations compare well with each 
other after scaling the MC times (expressed in number of sweeps) by a single multiplicative 
factor of 100. This shows that the relaxation dynamics does not depend on the details of 
the microscopic dynamics, 0s (t) exhibits an Arrhenius behavior with an activation energy 
of 0.28. 
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One of the most useful quantities for the study of dynamics is the intermediate scattering 
function F(k,t) defined as^ 

N N 

F ^ = ^ED ex P[-ik- -*i(0))]>- (11) 

1=1 3=1 

F(k,t) provides information on the collective dynamics of the system. The self intermedi- 
ate scattering function, F s (k,t), which reveals information about single particle motion, is 
obtained by restricting the double summation above to I = j: 

i N 

F s (k,t) = -^(exp[-ik- (r,(t) -r,(0))]). (12) 
j'=i 

Figure [10] shows F(k,t) vs. time t for different wave-vectors k at T — 0.04 and density 
p = 0.06. From this graph we recognize that F(k,t) shows three regimes: At very short 
times the function is quadratic in t, since we have a Newtonian dynamics. At long times 
the function shows at small k an oscillatory behavior which is related to the usual acoustic 
sound modes. The most remarkable feature is seen at intermediate times where we find that 
F(k,t) shows a compressed exponential decay, F(k,t) = exp( — (t/r)^), with an exponent 
f3 w 3/2 for all k- values shown. (For k values that are even larger than the ones shown in 
Fig. dm one finds that F(k,t) decays ballistically, i.e. ~ 2, see Ref.— .) Also included in 




FIG. 9: The bond lifetime tb vs. inverse T at density 0.06 and 0.10 from MD simulations and 
from MC simulations for density 0.06. 
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the graph are fits with a Kohlrausch- Williams- Watts (KWW) function with an exponent of 
3/2. Such compressed exponential relaxation, i.e. > 1, has been observed in experiments 
on non-equilibrium colloidal gels and it has been argued that this compression is related 
to the stress inhomogeneities caused by the shrinking of the gels^^£. In the present case, 
the system is in equilibrium and the compressed exponential behavior is due to the floppy 
motion of the non-restructuring network present in the gel-former—, as we will show below. 

It is also of interest to note that even at this low temperature there is no sign of the 
two-step relaxation observed in dense glass-forming liquids (see, e.g., Ref.— ). Thus we have 
here an example of a system whose relaxation dynamics is glassy (non-exponential, strong 
T— dependence, complex k— dependence,...) but with correlation functions that do not show 
the caging regime at intermediate times (at least in the k and T— range accessible in this 
simulation). 

In order to get a better understanding of the origin of the compressed exponential it 
is useful to compare the time dependence of F(k, t) with the one of the bond correlation 
function (j>B (t). This is done in Fig. [UJ, where we show F(k,t) for the smallest wave- vector 
as well as the bond correlation function for T = 0.04, p = 0.06. Also included is the bond 
correlation function defined only for particles that initially have three-fold coordination. It 
is seen that F(k, t) decays with a compressed exponential form on a time scale that is much 
shorter than the bond life time for all particles, and also in which a substantial fraction of 
bonds associated with initially three-coordinated particles are still intact. This observation 
suggests that the compressed exponential relaxation seen is due to the dynamics of the gel 
network on time scales when bond breaking is not very relevant. 

To confirm this, we have performed constrained MD simulations, i.e. a dynamics in which 
bonds are prevented from breaking or forming^. Figure [T2l shows F(k,t) with and without 
the presence of this constraint potential. The two F(k, t) curves are essentially identical in 
the time regime where compressed exponential relaxation is seen, but deviate from each other 
at later times, with the constrained MD curve saturating to a finite value. The data shown 
in Fig. [12] clearly demonstrates that the origin of the compressed exponential relaxation is 
in the vibrational dynamics of the gel network. 

Figure [TBI shows the time dependence of the self intermediate scattering function F s (k,t) 
at density 0.06 and temperature T = 0.04 for different wave-vectors. At long times, the 
decay is close to exponential for small wave vectors and faster than exponential at large 
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FIG. 10: The collective intermediate scattering function F(k,t) at p = 0.06 and T = 0.04 for 
different wave- vectors k. The KWW fit to F{k,t) at intermediate times with an exponent (3 = 3/2 
is also included (dashed lines). 

wave vectors^. A close inspection of Fig. [13] reveals that the stretching exponent varies 
non-monotonically, and is smaller at intermediate wave vectors. This behavior holds at 
other low temperatures as well^. A comparison with the behavior of F(k, t) shows that 
these two functions indeed reveal very different aspects of the dynamics of the gel former. 

In Fig. [14] we show the relaxation time r s that have been calculated from the area 
under F s (k,t) vs. time for different wave vectors at densities 0.06 (open symbols) and 
0.22 (filled symbols). (In order to compensate the trivial 1/T a5 dependence at large k, 
we multiply the data by T 0,5 .) At high wave-vectors the linear behavior of r s k 2 indicates 
the ballistic character of the dynamics, consistent with the gaussian behavior of F s (k,t) 
shown in Fig. [13] and the value for the KWW exponent which is around 2.0 (see Ref.— ). For 
intermediate temperatures, we observe a crossover to a regime where r s k 2 is roughly constant 
at small wave vectors, seemingly indicating a diffusive dynamics. However, the F s (k,t) 
curves shown in Fig. [13] reveal that the dynamics is more complex than a simple diffusive 
behavior in these cases. Note that the location in k at which this crossover is observed 
depends strongly on temperature and density, which indicates the complex dependence of 
the relaxation dynamics on p and T. This conclusion is also supported by the observation 
that at the lowest temperature investigated the mentioned plateau at low k disappears (see 
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Fig. IT4"|) since the stiffening of the network pushes the hydro dynamic regime to even lower 
values of k, in agreement with previous results^. 

The relaxation times (r and t s ) and tb at T = 0.04 and wave- vector k = 2tt x 2/L, 
i.e. the second smallest wave-vector compatible with the size L of the simulation box, is 
shown in Fig. [T5J The relaxation times t s , t have been obtained from the area under F s (k, t), 
F(k, t) vs. time curves. We recognize that tb as well as r decrease with increasing density. 
The reason for this is that with increasing p the number of short lived three and higher 
coordinated particles increases, with the consequences that: (i) tb of the system decreases 
and (ii) the network structure can reconstruct more easily, in turn leading to a faster decay 
of F(k,t), and a decrease of r. 

On the other hand, r s is nearly independent of density. This is a consequence of the 
fact that the MSD is roughly the same for all densities in the diffusive regime (see Fig. [6]). 
However, at larger wave-vectors, t s shows a more significant dependence on density (data 
not shown)^, consistent with significant changes in the MSD near the second shoulder. 

Since in real experiments the relaxation dynamics of the gels is given by a Brownian 
dynamics instead of the Newtonian dynamics studied here, it is of interest to investigate 
which properties of the relaxation dynamics depend on the microscopic dynamics. In the 




t 

FIG. 11: The collective intermediate scattering function F(k,t) and the bond correlation function 
<f>B(t) at density 0.06 and T = 0.04. Also shown is the bond correlation function for particles which 
are initially three- fold coordinated. 
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FIG. 12: The collective intermediate scattering function F(k,t) for different wave-vectors k at 
density p = 0.06 and T = 0.04 from MD (solid lines) and constrained MD simulation (dashed 
lines). 

final part of this paper we therefore compare F(k,t) and F s (k,t) obtained from MD with 
those obtained from MC. 

In Fig. [16] we show the coherent intermediate scattering function F(k,t) vs. time for 
different T's and k = 0.31 obtained from MC as well as MD simulations. In order to allow 




t 



FIG. 13: The self intermediate scattering function F s (k,t) from MD simulations for p = 0.06 and 
T = 0.04 for different wave-vectors k. 
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FIG. 14: The dependence of t s multiplied by k 2 , where r s has been obtained from the area of 
F s (k,t), at different T at p = 0.22 (filled symbols) and p = 0.06 (open symbols). 

to show the curves for different temperatures on the same graph, we plot the correlators as 
a function of t/t e , where t e is the time at which the F(k, t) reaches a value of 1/e. From the 
graph we see that the relaxation behavior of F(k,t) from the MC and MD simulations are 
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FIG. 15: The density dependence of the bond life time tb and relaxation times r and r s obtained 
from the area of the coherent and incoherent intermediate scattering functions for T = 0.04 and 
wave-vector k = ^ x 2L, where L is the size of the simulation box. 
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quite different: Independent of k^-, the F(k,t) from MC shows an exponential relaxation 
at the highest temperature shown (T = 0.10), and becomes progressively more stretched as 
temperature decreases. On the other hand, MD simulations exhibit compressed exponential 
for all T, although a slower decay is apparent at lower T due to the presence of acoustic 
modes. Thus we conclude that the MC dynamics suppresses the intermediate time com- 
pressed exponential behavior and allows one to see the stretched exponential decay at long 
times (see also Ref.— ). 

Last but not least we show in Fig. [T7] the incoherent intermediate scattering function 
F s (k,t) vs. t/t sc for different T and k obtained from MC as well as MD simulations. The 
F s (k,t) curves from MC are nearly exponential at the lower k values and higher T, showing 
stretching at the lowest T for all k values, and for all T at k = 1.23^. The degree of 
stretching increases with decrease of T and an increase of k. In contrast to this, F s (k,t) 
from MD shows a compressed exponential behavior for all T at k = 1.22^, but at smaller k 
the behavior becomes stretched at the two lowest temperatures. 




t/x e 

FIG. 16: The collective intermediate scattering function F(k, t) vs. time, for density = 0.06, in 
log-linear scale as obtained from MD (dashed lines) and MC (solid lines) for k = 0.31. 

Figure [18] shows the T— dependence of r obtained from the area under F(k, t) from MD 
and MC simulations for wave- vector k = 0.31, 0.61 and 1.22 at density 0.06. The data in Fig. 
[I8]have been plotted by rescaling the MC relaxation times to be the same at k = 0.31 for the 
highest T. We see that at T = 0.1 the A;— dependence of r is much stronger for the MD than 
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FIG. 17: The self intermediate scattering function F s (k, t) vs. time, for density = 0.06, in log-linear 
scale as obtained from MD (dashed lines) and MC (solid lines) for k = 0.31. 

for the MC. Furthermore we recognize that the T— dependence of r is stronger in the case of 
the MD than for the MC. Thus, the long time behavior as revealed by MC is substantially 
different from the self and collective intermediate scattering functions obtained from MD. 
This result also implies that from the point of view of simulations it is more efficient to 
equilibrate the system using MC. 



V. SUMMARY AND CONCLUSIONS 



In this paper we have investigated the static and dynamic properties of a recently intro- 
duced model for gel-forming systems. The gel-forming ability of this model is related to the 
fact that its local three-body potential avoids the formation of dense local structures and 
instead favors an open network structure. Due to this feature, the liquid-gas phase sepa- 
ration is pushed to low densities and temperatures, thus opening at intermediate densities 
a temperature range in which the relaxation dynamics is very slow, i.e. the system is a 
gel in equilibrium. Although this mechanism for the formation of the gel is similar to the 
approaches proposed earlier-^ 1 ^, the simplicity of the chosen interaction potential makes the 
present model very attractive for simulations. This advantage is enhanced even further by 
the use of an efficient method, presented in the Appendix, to evaluate the mentioned three- 
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FIG. 18: The T— dependence of the relaxation times obtained from the area under the coherent 
intermediate scattering function in MC (filled symbols) and MD (open symbols) simulations at 
density p = 0.06 for different wave-vectors. The MC data are scaled by a factor 2000 chosen such 
that the MC relaxation time at k = 0.31 coincides with the corresponding MD time at the highest 
temperature. 

body potential. Making use of this computational gain we have been able to characterize in 
detail the properties of the system, in equilibrium, in a significant range of temperature and 
density. 

We find that with decreasing temperature the particles assemble in the form of long one- 
dimensional chains that are connected to each other at random points, i.e. the system does 
indeed form a very heterogeneous network in which the length of the bridging chains follows 
an exponential distribution. The static structure factor of this spanning open network shows 
a strong increase at small wave- vectors for temperatures that are close to, but not quite at, 
the percolation line. The reason for this shift, which depends on density, is the fact that at 
large scales S(k) is not only given by the open structure of the percolating network but that 
there are also contributions from the form-factor of the chains. 

The bond correlation function 0s (t) shows at moderately low temperatures and high 
densities a non-exponential time dependence, which is related to the complex relaxation 
dynamics on the length scale of two particles, in analogy to the relaxation dynamics of 
dense glass-forming liquids. In contrast to this, 0# (t) shows at very low temperatures and 
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densities an exponential t— dependence, indicating a very simple bond breaking mechanism 
with an activation energy that is only a very weak function of T. We note that this result is 
independent of the used microscopic dynamics (molecular dynamics or Monte Carlo), which 
shows that this dynamical property is closely related to the structure of the system, as it 
has been found in other glass-forming systems^ - — . 

In contrast to this, the time correlation functions F(k,t) and F s (k,t) show a significant 
dependence on the microscopic dynamics. For F(k,t) as obtained from the MD we find a 
compressed exponential for length scales that are smaller or comparable to the one of the size 
of the particles, and acoustic modes for larger scales. On the other hand, the MC dynamics 
shows stretched exponentials for all length scales. The same trend is seen for F s (k, t), except 
that there are (of course) no acoustic modes in the MD. Thus we conclude that a dissipative 
dynamics, relevant for the experimental systems, does not show the compressed exponential 
relaxation observed in experiments in which the sample was aging. In experiments, therefore, 
the out-of-equilibrium nature of the samples must play an important role for the compressed 
exponential relaxation. Investigating the details of how the out-of-equilibrium dynamics 
leads to the experimentally observed compressed exponential relaxation remains an open 
problem that should be investigated in the future. 
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Appendix: Efficient method for force calculation of the SW potential 

Here we describe an efficient method for the calculation of energy and forces for the 
Stillinger- Weber potential, which involves only double sums to calculate three-body inter- 
action energies and forces. We follow the approach by Makhov and Lewis^ for energy 
calculations (see also^ 1 ^) but extend it to the calculation of forces. 



A. Energy Calculation 

From Eqs. ((!])- (jlj) it follows that the Stillinger- Weber potential can be written as 

u sw = J2Yl + Yl ^( r *i' r ^ ( 13 ) 

i j>i i jjti k^i,j 

with distances expressed in units of a and <fi denotes 



r ik ) = eX exp I 1 ) (cos 9 jik + a) 2 

\rij - a — a/ 

^H(a-r ij )H{a-r ik ) (14) 



We define A = eA, and the quantities 



9ij 



exp 




Tij < a 



nj > a 



(15) 



and 



A ^2 -gij g ik (cos9 jik + a) 2 

A ^2 [-g^ git® 2 + agij g ik cos 9 



(16) 



j=£i,k=£i 

+ gik cos 2 9 jik 



(17) 



With these definitions the three-body potential term in Eq. ( 1131) can be written as a sum 



over Uf 



25 



(18) 



where the correction Ufj' is due to the fact that Uf contains the terms j = k, see Eq. ( fl6l) . 
whereas the sum on the left hand side of Eq. f fl8|) does not include them. 
One easily finds that 



i.e. this is a two-body term. 
Defining thus 



1 - 



a 



(19) 



eff 
V- ■ 



Uf 



(20) 



where Uf, 



2Uf/, the total potential energy becomes 



(21) 



i j>i 



The factor 1/2 in Eq. (fl9l) is to account for the double counting of pair distances, which 
is not present in the first term of Eq. (]2ip and is, therefore, omitted. In the following we 
will show that the term Uf can be evaluated by a single loop over the neighbors of particle 
i, thus allowing us to evaluate the energy of the system without evaluating a triple sum. For 
this we define the quantity hi as 



and find 



hi 



(22) 



^ 9ij 9ik 



(23) 



which is proportional to the first term on the right hand side of Eq. (TTTI) . 
Similarly we define 



(24) 
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where r^- = ry/r^ to obtain 



Finally, we define T* 



9ij 9ik cos9 jik . 



^2 9v 9ik(iij-hk) 



(25) 
(26) 



T i = ^2g ij (r lj ®r lj ) 
where ® denotes the outer product, and obtain 



(27) 



Tr[T t 2 ] = ^ Tr[^-(fy (8) r^)^fe(fi fe ® r ife )] (28) 
= ^2 9ij9ik(iij-rik) 2 

= ^2 9ij9ik cos 2 9 jik . (29) 

Thus we find that Uf is a sum of three terms, each of which involves only a single sum over 
the neighbors of particle i: 

Uf = ^aX + Xa\s 2 \ + ^Tr[T 2 } (30) 

This expression thus allows us to calculate efficiently the total potential energy of the 
system using a double sum. 



B. Force Calculation 



The force acting on the i particle is given by 



dv 2 (rij) 
dra 



f)U c 



(31) 
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From Newton's third law of motion, we have 



Thus, 



., . ul ij ., . 



where derivatives of trivial two-body terms are given as 

-eA 



dv 2 (r i:j ) 



dri 



4Br' 5 H ^ 

13 (r y - a) 2 



x exp 



and 



2A(l + a) 2 ^^ 



The non-trivial part of the calculation is the calculation of VjUf 



Each of the terms are evaluated below. 



1. Calculation of Vj/i 2 : 



VA 2 



2h l V j h l 



2h i V j y^ j g ij 



2hi 



dgjj „ 
dri 



2 fei %j 
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2. Calculation of Vj | | : 

We have Vj|sf| = Vj(sj.Sj). The gradient of the dot product of two vectors X and Y 
is given by, 



V(X.Y) = (X.V)Y + (Y.V)X + X x (V x Y) 

+ Y x (V x X). (41) 

Therefore, Vj(sj.Sj) = 2(Vsi)sj. The gradient of vector s« is a tensor of rank two i.e. a 
matrix and is given as, 



VjSi = Vjigijiij) (42) 
- (r i:j <g> f + gijVj (r i:j ) (43) 



dr i3 



13 (iij ® Tij) + 9ji [I - Ti, ® fy] , (44) 



where we have used the fact that, 



VjTij = — (I - <g> fy). (45) 



Hence, 



= (46) 



Now, the gradient of |sj| is given by, 



2 



Vjls.l 2 = 2(V jSi ) Si (47) 
2 (dg ll _9 1 i] {hj ^ g . 



2^s, (48) 

T ■ ■ 



where we have used the fact that, 
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and, (iij ® r ij -)s l 



Finally, the gradient of s$ is written as 



2 ( r ij'- S iJ r ij- 



r 2 
ij 



Calculation of V^TrfT 2 
We use the fact that 
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dx n d(xA 
dx dx 



n-V 



A=x 



Therefore, for a constant matrix A, 



VjTr[T 2 ] = 2V i Tr(T i A 



A=Ti 

2V j Tr[g ij (r ij ® r -A)] 



One easily finds that, 



Trlgijiiij <g> r A)] 



9ij r ij 



Therefore, the gradient of trace of Tf matrix is given by 
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9ij^ij • Ar^ 



r ■ 

2@9ij ( T ij- r ^i T ij 



^ij ■ Ar j j 



A=T, 



1 U 



2 &9ij f ^-ij • i^-i] 



d rij 



1 ' 



1 y 



-) (rjj. Arjj)f j 



~ 2 9ij{~ 3 ; 



A=Ti 



A=T, 



Using Eq. (|4T|) . it is straightforward to show that 



V ; (r, r Ar,,) = 2Ar ?; 



because curl of and Ar^ vanish, i.e., 



V X Yij = 



V x Ar^ = (since A is a symmetric matrix). 



Therefore, finally the gradient of trace of Tf matrix becomes 



VjTr 



2 ^9ij ( r ij- r ^i r ij 



-Ag, 



ij 



r 



ij 



-2^(2T 4 r y ) 
T r 

+4q.. 4 ^ 



r, ; .T,r,j 



31 



Substituting Eqs. (SO}, ([52]), and §^ into Eq. ([36]), we get 



where, 



9ij\ ( ^i-^-ij \ ^9 



i j -\- Si> -\- 

ij 

c ij r ij + 2Aa — - Sj + 2\gij — 



dgij hA ~ ( dgij gij 



= [ ) +2\a ( - ) ( ^ 



A 



&9ij 2gij \ /Yij.TiYij 

dnj rij J V rfj 



The total force on a particle is obtained from Eq. ( 133]) as: 



£ (ViUf - v s uf) 



^2\g 

where 



(T i + Tj) 
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) 



r ii .(T i + T j ) r 



(70) 



This expression thus allows us to calculate efficiently the force on a particle, requiring 
only double sums. 
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Supplementary Material 



SI. CHOICE OF INTERACTION PARAMETERS a AND A 

We describe the procedure by which the interaction parameters a and A have been chosen. 
Figure [SU shows the LG phase-coexistence curves, obtained from Gibbs Ensemble Monte 
Carlo (GEMC) simulations^, for various combinations of A and a. The GEMC simulations 
have been performed with 2000 particles, with maximum step size of 0.085. The values 
A = 21 and a = 1/3 represent silicon and the corresponding phase-coexistence curve is the 
same obtained by Honda et al.—. When A is increased to 25 at constant a (= 1/3), we find 
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FIG. SI: The phase coexistence curves obtained from Gibbs Ensemble Monte Carlo (GEMC) 
simulations in T — p plane for (a) A=21, o=l/3, (b) A=25, o=l/3, (c) A=25, o=l/2, (d) A=10, 
a=1.0 and (e) A=10, a=l.l. By increasing the value of A or a, the phase coexistence region shifts 
to low T and p. The percolation line is also shown. The shaded region shows the estimated phase 
coexistence region for the present model for A=10, o=1.49 based on MD simulations (see text for 
details). The stars indicate the temperatures (at constant density) at which the structure factor 
at the smallest wave vector shows a maximum (see Fig. 2 of the main paper). The cross indicates 
the density and temperature at which the g(r) shown in Fig. [S2]is calculated. 
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that the LG phase-coexistence curve moves to lower temperatures. Further, if a is increased 
to 1/2 at fixed A (=25), the LG phase-coexistence curve shifts to even lower T's. If we keep 
A fixed at 10 and increase a, the LG phase coexistence curve also moves to low T and p 
values. Thus increasing either A or a results in diminishing the LG coexistence region. 

In addition to the effect on the LG coexistence, in order to ensure that the state of the 
system at intermediate densities is stable with respect to ordered structures (in our case, 
stackings of graphene-like sheets), we choose (among other possible choices) values A = 10 
and a = 1.49 (further details may be found in^). 

The LG phase coexistence curve for the values A = 10 and a = 1.49 can not be determined 
using GEMC simulations, because at low T bonds become too strong to swap particles 
between two sub volumes, which is required for GEMC simulations. However, from Fig. [SU 
it is apparent that for A = 10 and a = 1.49, the LG phase-coexistence curve will be further 
suppressed compared to A = 10 and a = 1.10. A possible way to deduce the location of 
the LG coexistence region is to look for phase separation in a constant volume simulation. 
The shaded area in Fig. [Si] shows the approximate region of LG phase coexistence for the 
present model from such MD simulations. However, the densities of both phases are very 
low, and therefore phase separation is not very easy to determine from the MD snapshots. 
Nevertheless, phase separation can be deduced by studying the radial distribution function, 
g(r), which, if the system phase separates, approaches to 1.0 from above instead of oscillating 
around 1.0 at large distances, the behavior found in homogeneous systems. This behavior 
is shown in Fig. IS2l 

For completeness, we show in Fig. [Si] also the percolation line. For a given temperature, 
the percolation density is determined as the density at which the percolation probability, 
estimated from considering several independent configurations, is 0.5 (i.e. 50% of the con- 
figurations have a spanning cluster of bonded particles). At the density p = 0.06, which we 
study in detail, the percolation transition occurs at T = 0.115. 

SII. STATIC PROPERTIES 

In order to get a first idea of the structure of the fluid, we show in Fig. [S3] snapshots from 
the MD simulations at density 0.06 and different temperatures. At a very high temperature, 
T = 5.0 (Fig. IS3h) the system is mainly composed of small sized clusters. Upon lowering the 
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FIG. S2: The radial distribution function, g(r), at density p = 0.009 for T = 0.035 (indicated 
by a cross in Fig. IS 1 P . The g(r) approaches its large-r limit from above, which indicates phase 
separation. 

temperature, particles form an increasing number of interconnections, forming interlinked 
linear chains that percolate if T is below 0.115. Figure IS3b shows a snapshot at T = 0.08 
manifesting large heterogeneity, and the presence of significant voids, due to the proximity 
to the percolation point. At still lower temperatures, particles form increasingly longer, and 
stiffer, chains, and the system displays a more homogeneous morphology (see Fig. IS3b for T = 
0.028). Therefore, at the lowest temperatures, the system consists mainly of linear chains 
of particles, with a small number of three- fold coordinated ( "anchor" ) particles connecting 
the linear chains. 

SIII. WAVE- VECTOR AND TEMPERATURE DEPENDENCE OF DYNAMICS 

In Fig. [S4] we show the coherent intermediate scattering function F(k,t) vs. time for 
different T and k values obtained from MC as well as MD simulations. In order to allow 
to show the curves for different temperatures on the same graph, we plot the correlators as 
a function of t/t e , where t e is the time at which the F(k,t) reaches a value of 1/e. From 
the graph we see that the relaxation behavior of F(k, t) from the MC and MD simulations 
are quite different: Independent of k, the F(k,t) from MC shows an exponential relaxation 
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FIG. S3: Snapshots of the system from MD simulations at density 0.06 for: (a) T = 5.0, (b) 
T = 0.08, and (c) T = 0.028. 

at the highest temperature shown (T = 0.10), and becomes progressively more stretched as 
temperature decreases. On the other hand, MD simulations exhibit compressed exponential 
for all T, although a slower decay is apparent at lower T due to the presence of acoustic 
modes. Thus we conclude that the MC dynamics suppresses the intermediate time com- 
pressed exponential behavior and allows one to see the stretched exponential decay at long 
times (see also Ref.-). 

Last but not least we show in Fig. [S5] the incoherent intermediate scattering function 
F s (k,t) vs. t/t se for different T and k obtained from MC as well as MD simulations. The 
F s (k, t) curves from MC are nearly exponential at the lower k values and higher T, showing 
stretching at the lowest T for all k values, and for all T at k — 1.22. The degree of 
stretching increases with decrease of T and an increase of k. In contrast to this, F s (k,t) 
from MD shows a compressed exponential behavior for all T at k — 1.22, but at smaller k 
the behavior becomes stretched at the two lowest temperatures (see main paper). 
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FIG. S4: The collective intermediate scattering function F(k, t) vs. time, for density = 0.06, in 
log-linear scale as obtained from MD (dashed lines) and MC (solid lines) for (a) k = 0.62 and (b) 
k = 1.22. 

SIV. CLUSTER SIZE DISTRIBUTION 

Fig. [S6]exhibits the cluster size distribution n s for different temperatures at density 0.06. 
For very high temperatures e.g., T = 5.0, the n s shows an exponential decay, consistent with 
the possibility of a a random aggregation process. When the temperature is lowered, the 
distribution develops a slower than exponential decay. However, even for temperatures near 
the percolation transition, we do not find a distinguishable power-law regime, in contrast to 
the findings for other gel forming systems for which a power-law dependence corresponding 
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FIG. S5: The self intermediate scattering function F s (k,t) vs. time, for density = 0.06, in log-linear 
scale as obtained from MD (dashed lines) and MC (solid lines) for (a) k = 0.62 and (b) k = 1.22. 

to random percolation has been observed^-. 
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FIG. S6: The cluster size distribution n s for different temperature at density 0.06. For very high 
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